#define c 137.037
#define pi 3.1416
#define nplim 2
#include <math.h>
#include <stdio.h>

void CalculaRo(double _xp[], int _np, int _npion, double _ro[], int _ndg,double _qp[],double _deltag, int _na, FILE *fp0) //na is a parameter of the sheet
{
	//int np,npion,ndg,
	//int i,j,na,l;
	//double xp[np],Ro[ndg],xg,qp[np],deltag,deltat,Tmax,qm;
	double xa,xb,qa,qb,qi,qil;
	double xg;
	double qtotalfinal, qtotalinicial;
	
	qtotalinicial=0;
    for(int j=0;j<_npion;j++)
        qtotalinicial+=_qp[j];
    //printf("***Carga total antes de calculaRo= %.4f\n", qtotalinicial);
	
	for (int i=0;i<_ndg;i++)
    {
		_ro[i]=0.;
    }
	
	for(int j=0;j<_npion;j++)
    {
		xg=_xp[j];
         
		//int i=int(xg/_deltag)+1;
        int i=int(xg/_deltag); //OJO!! MODIFICACION RESPECTO A FORTRAN, LE QUITO UNA UNIDAD PARA HACER EL GRID IGUAL
		//fprintf(fp0, "%d, %e, %e, %e\n", i, xg, xg/_deltag, (xg/_deltag)+1);
		//printf("\nPosicion=%.3f\tPunto en la red (i)=%d\tCarga=%.4f\n", xg, i, _qp[j]);
		//MODIFICACION PARA C: AUMENTO EN UNA UNIDAD LOS INDICES QUE MULTIPLICAN A DELTAG
        qa=_qp[j]*((i+1)*_deltag-xg)/_deltag;
		qb=_qp[j]*(xg-(i)*_deltag)/_deltag;
		//printf("qa=%.4f\tb=%.4f\n", qa, qb);
		
		xa=(xg+(i)*_deltag)/2.;
		xb=(xg+(i+1)*_deltag)/2.;
		qi=qa*((i+1)*_deltag-xa)/_deltag+qb*((i+1)*_deltag-xb)/_deltag;
		
		for (int l=i-_na/2;l<=i+_na/2;l++)
		{
			if(l>=1 && (l<=_ndg))
			{	_ro[l]=_ro[l]+qi/_na/_deltag;
                //printf("qi:Asigno %.4f al punto %d\n", qi/_na/_deltag, l);
            }
		}
		qil=qa*(xa-(i)*_deltag)/_deltag+qb*(xb-(i)*_deltag)/_deltag;
		//printf("qi=%.4f\tqi1=%.4f\n", qi, qil);
		for(int l=i+1-_na/2;l<=i+1+_na/2;l++)
		{
			if( (l>=1) && (l<=_ndg))
			{	_ro[l]=_ro[l]+qil/_na/_deltag;
				//printf("qil:Asigno %.4f al punto %d\n", qil/_na/_deltag, l);
            }
		}
		
	}
	

    /*for (int i=0;i<_ndg;i++)
        fprintf(fp0, "[%d]%e\t",i,_ro[i]);
    fprintf(fp0, "\n");*/
    
    /*printf("***Carga total tras calculaRo= %.4f\n", qtotalfinal*_deltag);*/
} // Calculate Ro

void CalculaRo_history(double _xp[], int _np, int _npion, double _ro[], int _ndg,double _qp[],double _deltag, int _na, FILE *fp0,int _n) //na is a parameter of the sheet
{
	//int np,npion,ndg,
	//int i,j,na,l;
	//double xp[np],Ro[ndg],xg,qp[np],deltag,deltat,Tmax,qm;
	double xa,xb,qa,qb,qi,qil;
	double xg;
	double qtotalfinal, qtotalinicial;
	
	qtotalinicial=0;
    for(int j=0;j<_npion;j++)
        qtotalinicial+=_qp[j];
    //printf("***Carga total antes de calculaRo= %.4f\n", qtotalinicial);
	
	for (int i=0;i<_ndg;i++)
    {
		_ro[i]=0.;
    }
	
	for(int j=0;j<_npion;j++)
    {
		xg=_xp[j];
        
		int i=int(xg/_deltag)+1;
		//printf("\nPosicion=%.3f\tPunto en la red (i)=%d\tCarga=%.4f\n", xg, i, _qp[j]);
		qa=_qp[j]*(i*_deltag-xg)/_deltag;
		qb=_qp[j]*(xg-(i-1)*_deltag)/_deltag;
		//printf("qa=%.4f\tb=%.4f\n", qa, qb);
		
		xa=(xg+(i-1)*_deltag)/2.;
		xb=(xg+i*_deltag)/2.;
		qi=qa*(i*_deltag-xa)/_deltag+qb*(i*_deltag-xb)/_deltag;
		
		for (int l=i-_na/2;l<=i+_na/2;l++)
		{
			if(l>=1 && (l<=_ndg))
			{	_ro[l]=_ro[l]+qi/_na/_deltag;
                //printf("qi:Asigno %.4f al punto %d\n", qi/_na/_deltag, l);
            }
		}
		qil=qa*(xa-(i-1)*_deltag)/_deltag+qb*(xb-(i-1)*_deltag)/_deltag;
		//printf("qi=%.4f\tqi1=%.4f\n", qi, qil);
		for(int l=i+1-_na/2;l<=i+1+_na/2;l++)
		{
			if( (l>=1) && (l<=_ndg))
			{	_ro[l]=_ro[l]+qil/_na/_deltag;
				//printf("qil:Asigno %.4f al punto %d\n", qil/_na/_deltag, l);
            }
		}
		
	}
	
	qtotalfinal=0;
    for (int i=0;i<_ndg;i++)
    {
		//fprintf(fp0, "%d\t%d\t%.5f\n",_n, i, _ro[i]);
        qtotalfinal+=_ro[i];
    }
    //printf("***Carga total tras calculaRo= %.4f\n", qtotalfinal*_deltag);
} // Calculate Ro


void CalculaJ(int _label, double _xp[],double _vpx[], double _vpy[],double _vpz[],int _np,int _npion,double _Jx[],double _Jy[],double _Jz[],int _ndg,double _qp[],double _deltat,double _deltag, int _na, FILE *fp0, FILE *fp1)
{
//	int np,npion,ndg,i,j,na.l;
//	double xp[np],vpx[np],vpy[np],vpz[np];
//	double Jx[ndg],Jy[ndg],Jz[ndg],xg,qp[ndg],deltag,deltat,Tmax,gamma,qm;
	double xa,xb,qa,qb,qi,qil,q,mass;
	double xg;
	
	if(_label==1)
	{
		for(int i=0;i<_ndg;i++)
			_Jx[i]=0.;
		
		for(int j=0;j<_npion;j++)
		{
			xg=_xp[j];
			//int i=int(xg/_deltag)+1; MODIFICACION PARA ADAPTACION A C
			int i=int(xg/_deltag);
			double gamma=sqrt(1.+(_vpx[j]*_vpx[j]+_vpy[j]*_vpy[j]+_vpz[j]*_vpz[j])/c/c);
			
            //MODIFICACION PARA C: AUMENTO EN UNA UNIDAD LOS INDICES QUE MULTIPLICAN A DELTAG
            qa=_qp[j]*((i+1)*_deltag-xg)/_deltag;
			qb=_qp[j]*(xg-(i)*_deltag)/_deltag;
			
			xa=(xg+(i)*_deltag)/2.;
			xb=(xg+(i+1)*_deltag)/2.;
			qi=qa*((i+1)*_deltag-xa)/_deltag+qb*((i+1)*_deltag-xb)/_deltag;
			
			for(int l=i-_na/2;l<=i+_na/2;l++)
			{
				if((l>=1) && (l<=_ndg))
					_Jx[l]=_Jx[l]+qi/_na/_deltag*(_vpx[j])/gamma;
			}
			
			qil=qa*(xa-(i)*_deltag)/_deltag+qb*(xb-(i)*_deltag)/_deltag;
			
			for(int l=i+1-_na/2;l<=i+1+_na/2;l++)
			{
				if((l>=1) && (l<=_ndg))
					_Jx[l]=_Jx[l]+qil/_na/_deltag*(_vpx[j])/gamma;
			}
			//fprintf(fp1, "%d\t%e\t%e\t%e\n",j, _vpx[j], _vpy[j], _vpz[j]);
			
		}//end j
		
    /*for (int i=0;i<_ndg;i++)
        fprintf(fp0, "[%d]: %e\n",i,_Jx[i]);*/
        
    //fprintf(fp0, "\n");		
       
		
	}//End lable1=Jx
	
	
	
	
	if(_label==2)
	{
		for(int i=0;i<_ndg;i++)
			_Jy[i]=0.;
		
		for(int j=0;j<_npion;j++)
		{
			xg=_xp[j];
			//int i=int(xg/_deltag)+1; //MODIFICACION PARA ADAPTACION A C
			int i=int(xg/_deltag);
			double gamma=sqrt(1.+(_vpx[j]*_vpx[j]+_vpy[j]*_vpy[j]+_vpz[j]*_vpz[j])/c/c);
			//MODIFICACION PARA C: AUMENTO EN UNA UNIDAD LOS INDICES QUE MULTIPLICAN A DELTAG
            qa=_qp[j]*((i+1)*_deltag-xg)/_deltag;
			qb=_qp[j]*(xg-(i)*_deltag)/_deltag;
			
			xa=(xg+(i)*_deltag)/2.;
			xb=(xg+(i+1)*_deltag)/2.;
			qi=qa*((i+1)*_deltag-xa)/_deltag+qb*((i+1)*_deltag-xb)/_deltag;
			
			for(int l=i-_na/2;l<=i+_na/2;l++)
			{
				if((l>=1) && (l<=_ndg))
					_Jy[l]=_Jy[l]+qi/_na/_deltag*_vpy[j]/gamma;
			}
			
			qil=qa*(xa-(i)*_deltag)/_deltag+qb*(xb-(i)*_deltag)/_deltag;
			
			for(int l=i+1-_na/2;l<=i+1+_na/2;l++)
			{
				if((l>=1) && (l<=_ndg))
					_Jy[l]=_Jy[l]+qil/_na/_deltag*(_vpy[j])/gamma;
			}
			
		}//end j
		
	}//end label2=Jy
	
	
	if(_label==3)
	{
		for(int i=0;i<_ndg;i++)
			_Jz[i]=0.;
		
		for(int j=0;j<_npion;j++)
		{
			xg=_xp[j];
			//int i=int(xg/_deltag)+1; //MODIFICACION PARA ADAPTACION A C
			int i=int(xg/_deltag);
			double gamma=sqrt(1.+(_vpx[j]*_vpx[j]+_vpy[j]*_vpy[j]+_vpz[j]*_vpz[j])/c/c);
			//MODIFICACION PARA C: AUMENTO EN UNA UNIDAD LOS INDICES QUE MULTIPLICAN A DELTAG
            qa=_qp[j]*((i+1)*_deltag-xg)/_deltag;
			qb=_qp[j]*(xg-(i)*_deltag)/_deltag;
			
			xa=(xg+(i)*_deltag)/2.;
			xb=(xg+(i+1)*_deltag)/2.;
			qi=qa*((i+1)*_deltag-xa)/_deltag+qb*((i+1)*_deltag-xb)/_deltag;
			
			for(int l=i-_na/2;l<=i+_na/2;l++)
			{
				if((l>=0) && (l<=_ndg))
					_Jz[l]=_Jz[l]+qi/_na/_deltag*(_vpz[j])/gamma;
			}
			
			qil=qa*(xa-(i)*_deltag)/_deltag+qb*(xb-(i)*_deltag)/_deltag;
			
			for(int l=i+1-_na/2;l<=i+1+_na/2;l++)
			{
				if((l>=0) && (l<=_ndg))
					_Jz[l]=_Jz[l]+qil/_na/_deltag*(_vpz[j])/gamma;
			}
			
		}//end j
		
	}//end label==3
	
}//End CalculaJ


double E0x(double _xg, double _t)
{
	//E0x=0.01*sin(2.*w0*t)*sin(wo*t)*sin(w0*t);
	double ex=0.;
	return ex;
}

double E0y(double _xg,double _t,double _xgleft, double _wp,double _kp,double _E0,double _k0,double _w0)
{
	double xx=_xgleft;
	double ey;
	
	if((_xg-c*_t < xx) && (_kp*(_xg-xx)-_wp*_t >=-pi) )
		ey=_E0*sin(_kp*(_xg-xx)-_wp*_t)*sin(_kp*(_xg-xx)-_wp*_t)*sin(_k0*(_xg-xx)-_w0*_t);
	else
		ey=0.;
	return ey;
}

double E0z(double _xg, double _t)
{
	double ez=0.;
	return ez;
}

double B0x(double _xg,double _t)
{
	double bx=0.;
	return bx;
}

double B0y(double _xg,double _t)
{
	double by=0.;
	return by;
}

double B0z(double _xg,double _t,double _xgleft, double _wp,double _kp,double _E0,double _k0,double _w0)
{
	double bz=E0y(_xg, _t, _xgleft,  _wp, _kp, _E0, _k0, _w0);
	return bz;
}


void CalculaEx(double _Ex[], double _Phi[], double _DAx[], int _ndg, double _t,double _deltag,double _deltat, FILE *fp0)
{
	// -grad Phi -1/c dAx/dt
	for(int i=1;i<(_ndg-1);i++)
	{
		//_Ex[i]=-(_Phi[i+1]-_Phi[i-1])/2./_deltag-_DAx[i]/c/_deltat+E0x((i-1)*_deltag,_t);//MODIFICACION DE FORTRAN A C
		_Ex[i]=-(_Phi[i+1]-_Phi[i-1])/2./_deltag-_DAx[i]/c/_deltat+E0x((i)*_deltag,_t);
	}
	_Ex[0]=-(_Phi[1]-_Phi[0])/_deltag-_DAx[0]/c/_deltat+E0x(0.,_t);
	_Ex[_ndg-1]=-(_Phi[_ndg-1]-_Phi[_ndg-2])/_deltag-_DAx[_ndg-1]/c/_deltat+E0x((_ndg-1)*_deltag,_t);
	
	/*for(int j=0; j<_ndg; j++)
	         fprintf(fp0, "%d\t%e\n", j, _Ex[j]);*/
}


void CalculaEy(double _Ey[],double _DAy[],int _ndg,double _t,double _deltag,double _deltat, double _xg,double _xgleft, double _wp,double _kp,double _E0,double _k0,double _w0)
{
	// -1/c dAx/dt
	for(int i=0;i<_ndg;i++)
	{
		//_Ey[i]=-_DAy[i]/c/_deltat+ E0y((i-1)*_deltag,_t,_xgleft, _wp,_kp, _E0,_k0, _w0); MODIFICACION PARA C
		_Ey[i]=-_DAy[i]/c/_deltat+ E0y((i)*_deltag,_t,_xgleft, _wp,_kp, _E0,_k0, _w0);
		//E0y( (i-1)*_deltag ,_t );
	}
}

void CalculaEz(double _Ez[],double _DAz[],int _ndg, double _t,double _deltag,double _deltat)
{
	// -1/c dAz/dt
	for(int i=0;i<_ndg;i++)
	{
		//_Ez[i]=-_DAz[i]/c/_deltat+E0z((i-1)*_deltag,_t); //MODIFICACION PARA C
		_Ez[i]=-_DAz[i]/c/_deltat+E0z((i)*_deltag,_t);
	}
}

void CalculaBy(double _By[],double _Az[], int _ndg, double _t,double _deltag,double _deltat)
{
	//-dAz/dx
	for(int i=1;i<(_ndg-1);i++)
	{
		//_By[i]=-(_Az[i+1]-_Az[i-1])/2./_deltag+B0y((i-1)*_deltag,_t); MODIFICACION PARA C
		_By[i]=-(_Az[i+1]-_Az[i-1])/2./_deltag+B0y((i)*_deltag,_t);
	}
	_By[0]=-(_Az[1]-_Az[0])/_deltag+B0y(0.,_t);
	//_By[_ndg-1]=-(_Az[_ndg-1]-_Az[_ndg-2])/_deltag+B0y((_ndg-1)*_deltag,_t); //MODIFICADO DE FORTRAN A C lo puse a 2, pero lo vuelvo a 1
	_By[_ndg-1]=-(_Az[_ndg-1]-_Az[_ndg-2])/_deltag+B0y((_ndg-1)*_deltag,_t); 
}

void CalculaBz(double _Bz[],double _Ay[],int _ndg,double _t,double _deltag,double _deltat,double _xg,double _xgleft, double _wp,double _kp,double _E0,double _k0,double _w0)
{
	//dAy/dx
	for(int i=1;i<(_ndg-1);i++)
	{
		//_Bz[i]=-(_Ay[i+1]-_Ay[i-1])/2./_deltag+B0z((i-1)*_deltag, _t, _xgleft,  _wp, _kp, _E0, _k0, _w0); MODIFICACION PARA C
        _Bz[i]=(_Ay[i+1]-_Ay[i-1])/2./_deltag+B0z((i)*_deltag, _t, _xgleft,  _wp, _kp, _E0, _k0, _w0);
   	}
	_Bz[0]=(_Ay[1]-_Ay[0])/_deltag+B0z(0.,_t, _xgleft,  _wp, _kp, _E0, _k0, _w0);
	//_Bz[_ndg-1]=-(_Ay[_ndg-1]-_Ay[_ndg-2])/_deltag+B0z((_ndg-1)*_deltag, _t, _xgleft,  _wp, _kp, _E0, _k0, _w0);//((_ndg-1),_t); MODIFICADO DE FORTRAN A C lo puse a 2 pero lo vuelvo a 1
	_Bz[_ndg-1]=(_Ay[_ndg-1]-_Ay[_ndg-2])/_deltag+B0z((_ndg-1)*_deltag, _t, _xgleft,  _wp, _kp, _E0, _k0, _w0);
	
}


void interpola(double _x,double _Ex[],double _Ey[],double _Ez[],double _By[],double _Bz[],double &_Epx,double &_Epy,double &_Epz,double &_Bpx,double &_Bpy,double &_Bpz, double _t,int _ndg,double _deltag, int _na, FILE *fp0)
{
	double xx=_x;
	//int i=int(xx/_deltag)+1;//MODIFICACION DE C
	int i=(int)round(xx/_deltag);
	_Epx=0.;
	_Epy=0.;
	_Epz=0.;
	_Bpx=0.;
	_Bpy=0.;
	_Bpz=0.;
	//if(i<309) printf("xx=%f, deltag=%f, xx/deltag=%f, (xx/deltag)+1=%f, int(xx/deltag)+1=%d\n", xx, _deltag, xx/_deltag, (xx/_deltag)+1, (int)(xx/_deltag)+1);
	//fprintf(fp0, "%d\t%e\txx/deltag=%e\t%e\t%e\t%e\t%e\t%e\t",i,_x, _x/_deltag, _Ex[i], _Ey[i], _Ez[i], _By[i], _Bz[i]);
	
	//printf("Int:");
	for(int l=i-_na/2;l<=i+_na/2;l++)
	{
		if((l>=0) && (l<_ndg) )
		{
			_Epx=_Epx+_Ex[l]/_na;
			_Epy=_Epy+_Ey[l]/_na;
			_Epz=_Epz+_Ez[l]/_na;
			
			_Bpx=_Bpx+B0x(xx,_t)/_na;
			_Bpy=_Bpy+_By[l]/_na;
			_Bpz=_Bpz+_Bz[l]/_na;
		}
	}
	//fprintf(fp0, "%e\t%e\t%e\t%e\t%e\t%e\n",_Epx, _Epy, _Epz, _Bpx, _Bpy, _Bpz);
	//printf("Epx=%.3f\tEpy=%.3f\tEpz=%.3f\tBpx=%.3f\tBpy=%.3f\tBpz=%.3f\t", *_Epx, *_Epy, *_Epz, *_Bpx, *_Bpy, *_Bpz);

}//end interpola

/*
void Ionizacion(int _ndg, int _nplim,double _deltap0[],double _p0[],double _Ex[],double _Ey[],double _Ez[],double _wrate,double _t,double _deltat,int _ileft,int _iright,double _coef,double _exponente)//,int _ppp[][]
{
	double wrate[_ndg];
	for(int i=_ileft;i<=_iright;i++)
	{
		double aux=sqrt(_Ex[i]*_Ex[i]+_Ey[i]*_Ey[i]+_Ez[i]*_Ez[i]);
		if(aux>3.333e-3)
		{
			wrate[i]=_coef/aux*exp(-_exponente/aux);
		}
		else
			wrate[i]=0.;
		
		_deltap0[i]=_deltat*wrate[i]*_p0[i];
		if(_deltap0[i]>_p0[i] )
			_deltap0[i]=_p0[i];
		//_ppp[i][1]=0;
	}
}


void CalculaP0(double _xp[],double _vpx[],double _vpy[],double _vpz[],int _np,int _npion, int _ndg,double _qp[],int _nplim, double _ppp[][nplim+1],double _p0[],double _rop[],double _deltap0[], int n, int _iseed,double _deltag,int _ileft, int _iright,double _cargaminima,int _nppp,int _npminimo,int _npsinp)
{
	int jal=1;//+int(_npion*rand(_iseed));
	for(int j=jal;j<_npion;j++)
	{
		int i=int(_xp[j]/_deltag)+1;
		if( (i*_deltag-_xp[j])< (_xp[j]-(i-1)*_deltag )    )
			i=i+1;
		if(_ppp[i][1] < _nplim )
		{
			_ppp[i][1]=_ppp[i][1]+1;
			_ppp[i][1]=j;//+_ppp[i][1]]=j;  //Check this line
		}
	}
	for(int j=1;j<jal-1;j++)
	{
		int i=int(_xp[j]/_deltag)+1;
		if( (i*_deltag-_xp[j]) < (_xp[j]-(i-1)*_deltag)  )
		{
			i=i+1;
		}
		if(_ppp[i][1]<_nplim )
		{
			_ppp[i][1]=_ppp[i][1]+1;
			_ppp[i][1]=j;// +_ppp[i,1]]=j;  fix this !!!
		}
	}
	
	for(int i=_ileft;i<_iright;i++)
	{
		if(_deltap0[i]> _cargaminima)
			_nppp=_ppp[i][1];
		if(_nppp!=0 && _npion>_npminimo)
		{
			_p0[i]=_p0[i]-_deltap0[i];
			_rop[i]=_rop[i]+_deltap0[i];
			double aux=_deltap0[i]*_deltag/_nppp;
			for(int j=1;j<_nppp;j++)
			{
				int jb=_ppp[i][1+j];
				double aux2=_qp[jb]-aux;
				if(aux2!=0)
				{
					_vpx[jb]=_vpx[jb]*_qp[jb]/aux2;
					_vpy[jb]=_vpy[jb]*_qp[jb]/aux2;
					_vpz[jb]=_vpz[jb]*_qp[jb]/aux2;
					_qp[jb]=aux2;
				}
			}//end for
		}//end if
		else if(_npion > _np)
		{
			_npion=_npion+1;
			_qp[_npion]=-_deltap0[i]*_deltag;
			_vpx[_npion]=0.;
			_vpy[_npion]=0.;
			_vpx[_npion]=0.;
			_xp[_npion]=_deltag*(i-1.);
			_p0[i]=_p0[i]-_deltap0[i];
			_rop[i]=_rop[i]+_deltap0[i];
		}
		else
			_npsinp=_npsinp+1;
	}
	
}


double ran1(int _idum,int _iff)
{
	int m=714025;
	int ia=1366;
	int ic=150889;
	double rm=1./m;
	double temp;
	int ir[97];
	int iy;
	
	if(_idum > 0 || _iff==0)
	{
		_iff=1;
		//_idum=mod(ic-_idum,m);
		for(int j=1;j<97;j++)
		{
			//_idum=mod(ia*_idum+ic,m);
			ir[j]=_idum;
		}
		//_idum=imod[ia*_idum+ic,m];
		iy=_idum;
	}
	int j=1+(97*iy)/m;
	//if(j>97 || j<1)
	//	break;
	
	iy=ir[j];
	double tempran1=iy*rm;
	//temp=imod(ia*_idum+ic,m);
	ir[j]=_idum;
	return tempran1;
}
*/

/*
void realft(double data,int n,int isignn)
{
	double pi=4.*atan(1.);
	double theta=pi/n;
	double cl=0.5;
	
	if(isign==1)
	{
		c2=-0.5;
		four1(data,n,1);
	}
	else
	{
		c2=0.5;
		theta=-theta;
	}
	double wpr=-2.*sin(0.5*theta)*sin(0.5*theta);
	double wpi=sin(theta);
	double wr=1.+wpr;
	double wi=wpi;
	n2p3=2*n+3;
	for(int i=2;i<n/2;i++)
	{
		int il=2*i-1;
		int i2=il+1;
		int i3=n2p3-i2;
		int i4=i3+1;
		double wrs=wr;
		double wis=wi;
		
		double hlr=c1*(data[i1]+data[i3]);
		double hli=cl*(data[i2]-data[i4]);
		double h2r=-c2*(data[i2]+data[i4]);
		double h2i=c2*(data[i1]-data[i3]);
		
		data[i1]=h1r+wrs*h2r-wis*h2i;
		data[i2]=h1i+wrs*h2i+wis*h2r;
		data[i3]=h1r-wrs*h2r+wis*h2i;
		data[i4]=-h1i+wrs*h2i+wis*h2r;
		double wtemp=wr;
		
		double wr=wr*wpr-wi*wpi*wr;
		double wi=wi*wpr+wtemp*wpi+wi;
	}
	
	if(isign==1)
	{
		h1r=data[1];
		data[1]=h1r+data[2];
		data[2]=h1r-data[2];
	}
	else
	{
		h1r=data[1];
		data[1]=c1*(h1r+data[2]);
		data[2]=c1*(h1r-data[2]);
		four1(data,n,-1);
	}
}


void four1(double data,int nn, int isign)
{
	double pi=4.*atan(1.);
	int n=2*nn;
	int j=1;
	
	for(int i=0;i<n;i++)
	{
			if(j>i)
			{
				tempr=data[j];
				tempi=data[j+1];
				
				data[j]=data[i];
				data[j+1]=data[i+1];
				data[i]=tempr;
				data[i+1]=tempi;
			}
			m=n/2;
			if(m>=2 && j > m)
			{
				j=j-m;
				m=m/2;
			}
			j=j+m;
			mmax=2;
			if(n>mmax)
			{
				istep=2*mmax;
				theta=2.*pi/(isign*mmax);
				wpr=-2.*sin(0.5*theta)*sin(0.5*theta);
				wpi=sin(theta);
				wr=1.;
				wi=0.;
				for(int m=0;m<mmax;m++2)
				{
					for(int i=m;i<n;i++step)
					{
						j=i+mmax;
						tempr=wr*data[j]-wi*data[j+1];
						tempi=wr*data[j+1]-wi*data[j];
						
						data[j]=data[i]-tempr;
						data[j+1]=data[i+1]-tempi;
						data[i]=data[i]+tempr;
						data[i+1]=data[i+1]+tempi;
					}
					wtemp=wr;
					wr=wr*wprwi*wpi+wr;
					wi=wi*wpr+wtemp*wpi+wi;
				}
				mmax=istep;
			}
	}
}

*/
